---
title: "DDR dataset"
author: "Sally Sharif"
date: "01/29/2026"
---

```{r loading packages and data, message=FALSE, warning=FALSE}
rm(list=ls()) 
# loading required packages
library(texreg)
library(dplyr)
library(ggplot2)
library(ggrepel)
library(scales)
library(stringr)
library(tidyr)
library(table1)
library(corrtable)
library(glmmTMB)
library(performance)
library(car)
library(rpart)
library(car)
library(rpart.plot)
library(randomForest)

# loading the data
obj_name <- load("DDR40.RData")
DDR40 <- get(obj_name)
nrow(DDR40) # 407
```

```{r recoding variables}
# clean + recode core variables (rearmament, budget measures, factor references, group types, geography)
DDR40 <- DDR40 %>%
  mutate(
    # rearmament during DDR: 1 if group OR any faction rearmed; 0 if neither; NA otherwise
    rearmament_during_ddr = case_when(
      rearmed == 1 | factions_rearmed == 1 ~ 1L,
      rearmed == 0 & factions_rearmed == 0 ~ 0L,
      TRUE ~ NA_integer_),
    # budget as numeric
    budget = as.numeric(budget),
    # budget per combatant (avoid dividing by zero)
    budget_per_combatant = budget / pmax(numbers, 1),
    # log-transform budget per combatant (add 1 so zeros are defined)
    log_budget_per_combatant = log(budget_per_combatant + 1),
    # set reference categories for key factor covariates
    establishment = relevel(factor(establishment), ref = "national initiative"),
    execution     = relevel(factor(execution),     ref = "National"),
    group_type_tri = relevel(factor(group_type_tri), ref = "Military"),
    # create a binary variable: state/non-state groups undergoing DDR
    group_type_binary = case_when(
      group_type_tri %in% c("Fragmented", "Integrated") ~ 1L,
      group_type_tri == "Military" ~ 0L,
      TRUE ~ NA_integer_),
    # binary non-state group type: Integrated=1, Fragmented=0, NA otherwise (incl. Military)
    nonstate_group_type = case_when(
      group_type_tri == "Integrated"  ~ 1L,
      group_type_tri == "Fragmented"  ~ 0L,
      TRUE                            ~ NA_integer_),
    # continent from sub-region labels
    continent = case_when(
      subregion %in% c("North Africa", "West Africa", "Central Africa", 
                       "East Africa", "Southern Africa",
                       "Great Lakes", "Horn of Africa") ~ "Africa",
      subregion %in% c("Central America", "Caribbean", "South America") ~ "America",
      subregion %in% c("South Asia", "Southeast Asia", "Central Asia", "West Asia",
                       "Middle East", "Eastern Asia", "Southwest Asia", 
                       "Pacific Islands") ~ "Asia",
      subregion %in% c("Europe", "Eastern Europe") ~ "Europe",
      TRUE ~ NA_character_)) %>%
  # add a scaled version of log budget per combatant (z-score)
  mutate(log_budget_per_combatant_scaled = as.numeric(scale(log_budget_per_combatant)))

# combined provision + implementation variables:
# rule: if the accord has no provision, set combined score to 0; otherwise keep the implementation score (1–3)
DDR40 <- DDR40 %>%
  mutate(
    amnesty_combined               = if_else(amnesty_provision == 0, 0L,
                                             as.integer(amnesty_implementation)),
    prisoner_release_combined      = if_else(prisoner_release_provision == 0, 0L,
                                             as.integer(prisoner_release_implementation)),
    truth_reconciliation_combined  = if_else(truth_reconciliation_provision == 0, 0L,
                                             as.integer(truth_reconciliation_implementation)),
    womens_rights_combined         = if_else(womens_rights_provision == 0, 0L,
                                             as.integer(womens_rights_implementation)),
    childrens_rights_combined      = if_else(childrens_rights_provision == 0, 0L,
                                             as.integer(childrens_rights_implementation)),
    boundary_demarcation_combined  = if_else(boundary_demarcation_provision == 0, 0L,
                                             as.integer(boundary_demarcation_implementation)),
    decentralization_combined      = if_else(decentralization_provision == 0, 0L,
                                             as.integer(decentralization_implementation)),
    dispute_resolution_combined    = if_else(dispute_resolution_provision == 0, 0L,
                                             as.integer(dispute_resolution_implementation)),
    constitutional_reform_combined = if_else(constitutional_reform_provision == 0, 0L,
                                             as.integer(constitutional_reform_implementation)),
    electoral_reform_combined      = if_else(electoral_reform_provision == 0, 0L,
                                             as.integer(electoral_reform_implementation)),
    executive_reform_combined      = if_else(executive_reform_provision == 0, 0L,
                                             as.integer(executive_reform_implementation)),
    judiciary_reform_combined      = if_else(judiciary_reform_provision == 0, 0L,
                                             as.integer(judiciary_reform_implementation)),
    military_reform_combined       = if_else(military_reform_provision == 0, 0L,
                                             as.integer(military_reform_implementation)),
    legislative_reform_combined    = if_else(legislative_reform_provision == 0, 0L,
                                             as.integer(legislative_reform_implementation)),
    development_combined           = if_else(development_provision == 0, 0L,
                                             as.integer(development_implementation)),
    natural_resources_combined     = if_else(natural_resources_provision == 0, 0L,
                                             as.integer(natural_resources_implementation)))
```

```{r descriptive statistics}
## Manuscript Table 2
# descriptive statistics 
table1(~ numbers + factor(multiple_groups) + 
         disarmament_implementation + demobilization_implementation + 
         reintegration_implementation + 
         cantonment + factor(establishment) + factor(execution) + 
         budget + factor(group_type_tri) + factor(unpko), data=DDR40) 

## Manuscript Table 3
# descriptive statistics of implementation for variables with the provision
variable_pairs <- list(
  c("amnesty_provision", "amnesty_implementation"),
  c("prisoner_release_provision", "prisoner_release_implementation"),
  c("boundary_demarcation_provision", "boundary_demarcation_implementation"),
  c("decentralization_provision", "decentralization_implementation"),
  c("detailed_timeline_provision", "detailed_timeline_implementation"),
  c("dispute_resolution_provision", "dispute_resolution_implementation"),
  c("development_provision", "development_implementation"),
  c("electoral_reform_provision", "electoral_reform_implementation"),
  c("constitutional_reform_provision", "constitutional_reform_implementation"),
  c("executive_reform_provision", "executive_reform_implementation"),
  c("judiciary_reform_provision", "judiciary_reform_implementation"),
  c("legislative_reform_provision", "legislative_reform_implementation"),
  c("military_reform_provision", "military_reform_implementation"),
  c("natural_resources_provision", "natural_resources_implementation"),
  c("truth_reconciliation_provision", "truth_reconciliation_implementation"),
  c("womens_rights_provision", "womens_rights_implementation"),
  c("childrens_rights_provision", "childrens_rights_implementation"))

# initialize an empty data frame to store the results
results_df <- data.frame(
  provision_var = character(),
  implementation_var = character(),
  mean = numeric(),
  sd = numeric(),
  stringsAsFactors = FALSE)

# loop through each pair and calculate the statistics
for (pair in variable_pairs) {
  provision_var <- pair[1]
  implementation_var <- pair[2]
  # check if both provision and implementation columns exist in DDR40
  if (all(c(provision_var, implementation_var) %in% colnames(DDR40))) {
    # filter rows where provision variable is 1
    filtered_df <- DDR40 %>% filter(!!sym(provision_var) == 1)
    # print the number of rows that match the filter for debugging
    cat("\n", provision_var, "and", implementation_var, " - Rows matching filter: ", nrow(filtered_df), "\n")
    # if there are matching rows, calculate the mean and SD
    if (nrow(filtered_df) > 0) {
      stats <- filtered_df %>%
        summarise(
          mean = mean(!!sym(implementation_var), na.rm = TRUE),
          sd = sd(!!sym(implementation_var), na.rm = TRUE))
      # append the results to the results_df
      results_df <- rbind(results_df, data.frame(
        provision_var = provision_var,
        implementation_var = implementation_var,
        mean = stats$mean,
        sd = stats$sd))}}}
# view the consolidated results
print(results_df)

## Manuscript Table 4
table1(~ factor(rearmed) + factor(factions_rearmed) + factor(rearmed_post_ddr), data=DDR40)

# Appendix Table 3
selected_data <- DDR40[, c("numbers", "multiple_groups", "disarmament_implementation", 
                           "demobilization_implementation", "reintegration_implementation", 
                           "cantonment", "budget", "rearmed", "factions_rearmed", 
                           "rearmed_post_ddr")]

correlation_matrix(selected_data, digits = 2, use = "lower", replace_diagonal = TRUE)
save_correlation_matrix(df=selected_data, filename = "corr.csv", digits=2, use="lower")

# Appendix table 4
# correlations
selected_data2 <- DDR40[, c("amnesty_implementation", "prisoner_release_implementation",
                           "truth_reconciliation_implementation", "womens_rights_implementation",
                           "childrens_rights_implementation", "development_implementation",
                           "dispute_resolution_implementation", "boundary_demarcation_implementation",
                           "natural_resources_implementation", "decentralization_implementation",
                           "constitutional_reform_implementation", "electoral_reform_implementation", 
                           "executive_reform_implementation", "judiciary_reform_implementation", 
                           "military_reform_implementation", "legislative_reform_implementation")]

correlation_matrix(selected_data2, digits = 2, use = "lower", replace_diagonal = TRUE)
save_correlation_matrix(df=selected_data2, filename = "corr.csv", digits=2, use="lower")
```

```{r descriptive graphs}
## Manuscript Figure 1
# summarize the data by continent
plot_df <- DDR40 %>%                       
  group_by(id) %>%                         
  summarise(start_year = min(year),
            combatants = first(numbers),
            continent  = first(continent),
            .groups = "drop") %>% 
  mutate(log_combatants = log10(combatants))

# find label candidates by continent
largest_ids <- plot_df %>% 
  group_by(continent) %>% slice_max(combatants, n = 1, with_ties = FALSE) %>% 
  pull(id)

recent_ids <- plot_df %>% 
  group_by(continent) %>% slice_max(start_year, n = 1, with_ties = FALSE) %>% 
  pull(id)

label_ids <- union(largest_ids, recent_ids)   
plot_df   <- plot_df %>% mutate(to_label = id %in% label_ids)

# draw the figure
p <- ggplot(plot_df,
       aes(x = start_year,
           y = log_combatants,
           colour = continent)) +
  geom_point(size = 3, alpha = .85) +
  geom_text_repel(data = subset(plot_df, to_label),
                  aes(label = id),
                  size = 3, box.padding = .25,
                  nudge_y = .15, family = "sans") +
  scale_colour_manual(values = c(
      "Africa"   = "#1b9e77",
      "America"  = "#d95f02",
      "Asia"     = "#7570b3",
      "Europe"   = "#e7298a"),
      name = NULL) +
  scale_y_continuous(name = "Target combatants for DDR (log)") +
  scale_x_continuous(name = "Program start year",
                     breaks = seq(1980, 2020, 10)) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom")

## save the figure at print‑quality resolution 
# ggsave(filename = "DDR40recent.png", plot = p, width = 7, height = 5, units    = "in")

## Manuscript Figure 2
# last year of each program
DDR_last <- DDR40 %>%
  group_by(id) %>% slice_max(year, n = 1, with_ties = FALSE) %>% ungroup()

# reshape & clean provision names 
clean_levels <- c(
  "Disarmament","Demobilization","Reintegration",
  "Amnesty","Prisoner Release","Dispute Resolution","Truth and Reconciliation",
  "Development","Natural Resources",
  "Boundary Demarcation","Decentralization","Detailed Timeline",
  "Constitutional Reform","Electoral Reform","Executive Reform",
  "Legislative Reform","Judiciary Reform","Military Reform",
  "Women's Rights","Children's Rights")

impl_vars <- paste0(
  c("amnesty","prisoner_release","boundary_demarcation",
    "decentralization","detailed_timeline","dispute_resolution",
    "disarmament","demobilization","reintegration",
    "development","electoral_reform","constitutional_reform",
    "executive_reform","judiciary_reform","legislative_reform",
    "military_reform","natural_resources","truth_reconciliation",
    "womens_rights","childrens_rights"),
  "_implementation")

# keep only final year of each programme
DDR_last <- DDR40 %>% group_by(id) %>% slice_max(year, n = 1) %>% ungroup()

impl_long <- DDR_last %>%
  select(all_of(impl_vars)) %>%
  mutate(rowID = row_number()) %>%                           
  pivot_longer(-rowID, names_to = "provision", values_to = "score") %>%
  filter(!is.na(score)) %>%                                  
  mutate(
    provision = provision %>% 
      str_remove("_implementation$") %>% 
      str_replace_all("_", " ") %>% 
      str_to_title() %>% 
      str_replace("Truth Reconciliation", "Truth and Reconciliation") %>% 
      str_replace("Womens Rights", "Women's Rights") %>% 
      str_replace("Childrens Rights", "Children's Rights"),
    provision = factor(provision, levels = clean_levels, ordered = TRUE))

# quick check: all names match?
setdiff(levels(impl_long$provision), clean_levels)   # should return character(0)

q <- ggplot(impl_long, aes(x = factor(score))) +
  geom_bar(fill = "#255C99", width = .85) +
  facet_wrap(~ provision, ncol = 4, scales = "free_y") +
  scale_x_discrete(drop = FALSE, limits = c("0","1","2","3")) +
  labs(x = "Implementation score (0–3)",
       y = "Number of programmes in final year") +
  theme_minimal(base_size = 9) +
  theme(panel.grid.minor = element_blank(),
        strip.text = element_text(size = 8))

# ggsave("implementation.png", q, width = 7, height = 5.2, dpi = 600)

## Manuscript Figure 3
# collapse the three binary outcomes into one categorical variable
DDR40_plot <- DDR40 %>% 
  mutate(establishment = factor(establishment,
           levels = c("peace agreement","national initiative","international initiative")),
         outcome = case_when(
           rearmed == 1                ~ "Group rearmed during DDR",
           factions_rearmed == 1       ~ "Faction rearmed during DDR",
           rearmed_post_ddr == 1       ~ "Rearmed after DDR",
           TRUE                        ~ "No rearmament"), 
        outcome = factor(outcome,
           levels = c("No rearmament","Group rearmed during DDR",
                      "Faction rearmed during DDR","Rearmed after DDR")))
# draw share‑of‑program‑years bars, one facet per establishment mode
r <- ggplot(DDR40_plot, aes(x = establishment, fill = outcome)) +
  geom_bar(position = "fill", width = .8, color = "grey10") +
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  scale_fill_brewer(palette = "Paired", name = NULL) +
  labs(x = NULL, y = "Share of program years") +
  coord_flip() +
  theme_minimal(base_size = 10) +
  theme(legend.position = "bottom")

# ggsave("DDR40_outcomes.png", r, width = 7.2, height = 3.5, dpi = 600)
```

```{r main models with program features}
## Manuscript Table 5
# pooled logit (Column 1)
pooled_logit1 <- glm(rearmament_during_ddr ~ I(numbers/1000) + multiple_groups + 
                       disarmament_implementation + demobilization_implementation + 
                       reintegration_implementation + log_budget_per_combatant + 
                       establishment + execution + unpko + log(cantonment+1), 
                     data = DDR40, family = binomial)
# check variance inflation
vif(pooled_logit1)

# pooled logit with binary group type: state vs. non-state (Column 2)
pooled_logit2 <- glm(rearmament_during_ddr ~ I(numbers/1000) + multiple_groups +
                       group_type_binary + 
                       disarmament_implementation + demobilization_implementation + 
                       reintegration_implementation + log_budget_per_combatant + 
                       establishment + execution + unpko + log(cantonment+1), 
                     data = DDR40, family = binomial)
# check variance inflation
vif(pooled_logit2)

# pooled logit with binary non-state group: integrated vs. fragmented (Column 3)
pooled_logit3 <- glm(rearmament_during_ddr ~ I(numbers/1000) + multiple_groups +
                       nonstate_group_type + 
                       disarmament_implementation + demobilization_implementation + 
                       reintegration_implementation + log_budget_per_combatant + 
                       establishment + execution + unpko + log(cantonment+1), 
                     data = DDR40, family = binomial)
# check variance inflation
vif(pooled_logit3)

# random effects logit with binary group type: state vs. non-state (Column 4)
random_effects1 <- glmmTMB(rearmament_during_ddr ~ I(numbers / 1000) + multiple_groups + 
                             group_type_binary +
                             disarmament_implementation + demobilization_implementation + 
                             reintegration_implementation + establishment + execution + 
                             unpko + log(cantonment+1) +
                             (1 | id), data = DDR40, family = binomial(link = "logit"))

# random effects logit with binary non-state group: integrated vs. fragmented (Column 5)
random_effects2 <- glmmTMB(rearmament_during_ddr  ~ I(numbers / 1000) + multiple_groups + 
                             nonstate_group_type +
                             disarmament_implementation + demobilization_implementation + 
                             reintegration_implementation + execution + 
                             unpko + log(cantonment + 1) +
                             (1 | id), data = DDR40, family = binomial(link = "logit"))

# fixed effects logit (Column 6)
fixed_effects <- glm(rearmament_during_ddr ~ I(numbers/1000) + multiple_groups + 
                 disarmament_implementation + demobilization_implementation + 
                 reintegration_implementation 
                 + as.factor(country), data = DDR40, family = binomial) 
# removed execution and cantonment and unpko because of fixed effects
AIC(fixed_effects)
logLik(fixed_effects)

# table
texreg(list(pooled_logit1, pooled_logit2, pooled_logit3, random_effects1,
            random_effects2, fixed_effects))

## testing whether cantonment and ex-combatant numbers need to be transformed
# ensure positive (add small constant if zeros exist)
DDR40 <- DDR40 %>%
  mutate(cantonment_pos = if_else(cantonment <= 0, min(cantonment[cantonment > 0], 
                                                       na.rm=TRUE)/2, cantonment))

# Box-Tidwell for numbers
model_bt_numbers <- glm(rearmament_during_ddr ~ numbers + I(numbers * log(numbers)) +
                         multiple_groups +  group_type_binary + 
                          disarmament_implementation + demobilization_implementation +
                          reintegration_implementation + log_budget_per_combatant +
                          establishment + execution + unpko,
                        family = binomial, data = DDR40)

summary(model_bt_numbers)
# no need for transformation of numbers: the logit-linearity assumption is not violated. 
# the Box–Tidwell interaction term I(numbers_pos * log(numbers_pos)) is tiny and far from significant.

# Box-Tidwell for cantonment
model_bt_canton <- glm(rearmament_during_ddr ~ cantonment_pos + I(cantonment_pos * log(cantonment_pos)) +
                         numbers + group_type_binary + multiple_groups +
                         disarmament_implementation + demobilization_implementation +
                         reintegration_implementation + log_budget_per_combatant +
                         establishment + execution + unpko,
                       family = binomial, data = DDR40)

summary(model_bt_canton)
# cantonment needs to be logged
```

```{r main models with peace provisions}
## Manuscript Table 6 and Appendix Table 5
# amnesty (Column 1)
random_effects3 <- glmmTMB(rearmament_during_ddr ~ I(numbers / 1000) + multiple_groups + 
                             execution + group_type_binary +
                             disarmament_implementation + demobilization_implementation + 
                             reintegration_implementation + as.factor(amnesty_combined) + 
                             establishment + log(cantonment+1) + (1 | id), 
                           data = DDR40, family = binomial(link = "logit"))

# prisoner release (Column 2)
random_effects4 <- glmmTMB(rearmament_during_ddr ~ I(numbers / 1000) + multiple_groups + 
                             execution + group_type_binary + 
                             disarmament_implementation + demobilization_implementation + 
                             reintegration_implementation +
                             as.factor(prisoner_release_combined) + 
                             establishment + log(cantonment+1) + (1 | id), 
                           data = DDR40, family = binomial(link = "logit"))

# dispute resolution (Column 3)
random_effects5 <- glmmTMB(rearmament_during_ddr ~ I(numbers / 1000) + multiple_groups + 
                             execution + group_type_binary + 
                             disarmament_implementation + demobilization_implementation + 
                             reintegration_implementation + 
                             as.factor(dispute_resolution_combined) + 
                             establishment + log(cantonment+1) + (1 | id), 
                           data = DDR40, family = binomial(link = "logit"))


# executive reform (Column 4)
random_effects6 <- glmmTMB(rearmament_during_ddr ~ I(numbers / 1000) + multiple_groups + 
                             execution + group_type_binary + 
                             disarmament_implementation + demobilization_implementation + 
                             reintegration_implementation +  
                             as.factor(executive_reform_combined) +
                             establishment + log(cantonment+1) + (1 | id), 
                           data = DDR40, family = binomial(link = "logit"))

# decentralization (Column 5)
random_effects7 <- glmmTMB(rearmament_during_ddr ~ I(numbers / 1000) + multiple_groups + 
                             execution + group_type_binary +
                             disarmament_implementation + demobilization_implementation + 
                             reintegration_implementation +
                             as.factor(decentralization_combined) + 
                             establishment + log(cantonment+1) + (1 | id), 
                           data = DDR40, family = binomial(link = "logit"))

# constitutional reform (Column 6)
random_effects8 <- glmmTMB(rearmament_during_ddr ~ I(numbers / 1000) + multiple_groups + 
                             execution + group_type_binary +
                             disarmament_implementation + demobilization_implementation + 
                             reintegration_implementation + 
                             as.factor(constitutional_reform_combined) + 
                             establishment + log(cantonment+1) + (1 | id), 
                           data = DDR40, family = binomial(link = "logit"))

# legislative reform (Column 7)
random_effects9 <- glmmTMB(rearmament_during_ddr ~ I(numbers / 1000) + multiple_groups + 
                             execution + group_type_binary +
                             disarmament_implementation + demobilization_implementation + 
                             reintegration_implementation +
                             as.factor(legislative_reform_combined) +
                             establishment + log(cantonment+1) + (1 | id), 
                           data = DDR40, family = binomial(link = "logit"))

# table
texreg(list(random_effects3, random_effects4, random_effects5, random_effects6,
            random_effects7, random_effects8, random_effects9))
```

```{r decision tree}
# single decision tree with program features
data_trees <- DDR40 %>% select(rearmament_during_ddr, numbers, multiple_groups, 
                               disarmament_implementation, demobilization_implementation, 
                               reintegration_implementation, cantonment, budget)
# Appendix Figure 1
mod.CART <- rpart(rearmament_during_ddr ~., data = data_trees, 
                  control = rpart.control(cp = 0.02))

# open a PNG device with high resolution
png("tree.png", width = 800, height = 800, res = 300)

rpart.plot(mod.CART, box.palette = "Grays", type = 2, 
           fallen.leaves = TRUE)
dev.off()
```

```{r random forest}
# selecting all implementation variables
data_forest <- DDR40 %>% select(rearmament_during_ddr, disarmament_implementation,
                                demobilization_implementation, reintegration_implementation, 
                                amnesty_implementation, prisoner_release_implementation,
                                boundary_demarcation_implementation, decentralization_implementation,
                                dispute_resolution_implementation, development_implementation,
                                electoral_reform_implementation, constitutional_reform_implementation,
                                executive_reform_implementation, judiciary_reform_implementation,
                                legislative_reform_implementation, military_reform_implementation,
                                natural_resources_implementation, truth_reconciliation_implementation,
                                womens_rights_implementation, childrens_rights_implementation)
# for replicability
set.seed(123)

# RF with peace provisions
rfmod <- randomForest(x = data_forest[-1], y = as.factor(data_forest$rearmament_during_ddr), 
                      importance = TRUE, ntree = 10000, maxnodes = 5, 
                      sampsize = 100, replace = FALSE, do.trace = FALSE, 
                      data = data1, forest = TRUE)   
importance(rfmod)
## Manuscript Figure 4
varImpPlot(rfmod)	
# high-res PNG, 300 dpi, showing only Mean Decrease Accuracy
png("varImp_MDA_300dpi.png", width = 6, height = 5.5, units = "in", res = 300)
dev.off()
```

```{r homicides}
# load the data with homicides
obj_name <- load("homicides.RData")
data_homicides <- get(obj_name)
nrow(data_homicides) # 407
## Appendix Table 6
# all groups (Column 1)
linear_mod1 <- lm(homicide_rate ~ I(numbers/1000) + multiple_groups + 
                    disarmament_implementation + demobilization_implementation +
                    reintegration_implementation + log_budget_per_combatant_scaled + 
                    establishment + execution + unpko + log(cantonment+1) +
                    factor(country) + factor(year), 
                  data = data_homicides)

AIC(linear_mod1)
# only non-state groups (Column 2)
linear_mod2 <- lm(homicide_rate ~ I(numbers/1000) + multiple_groups + 
                    group_type_binary + 
                    disarmament_implementation + demobilization_implementation +
                    reintegration_implementation + log_budget_per_combatant_scaled + 
                    establishment + execution + unpko + log(cantonment+1) +
                    factor(country) + factor(year), 
                  data = data_homicides)

AIC(linear_mod2)
# table
texreg(list(linear_mod1, linear_mod2))
```
